cd "<filepath>"

* Automatically generate new log file
local c_date = c(current_date)
local c_time = c(current_time)

local c_time_date = "`c_date'"+"_"+"`c_time'"
local time_str = subinstr("`c_time_date'", ":", "_", .)
local time_str = subinstr("`time_str'", " ", "_", .)
di "`time_str'"

log using "Replication_RI_`time_str'_NJ & VA 2015.txt", text 

*========================================
* VA & NJ 2015 GOTV Program
* -> Analysis <-
* -> Randomization Inference <-
*========================================

* Version of stata
version

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Set Directory
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

cd "<filepath>"


* Expand available matrix size
set matsize 10000

	
*========================================
* Set RI paramaters & postfile
*========================================

local va_overall_n "ri_va_eng_n ri_va_bil_n ri_va_diff_n" 
local va_overall_c "ri_va_eng_c ri_va_bil_c ri_va_diff_c" 

local nj_overall_n "ri_nj_eng_n ri_nj_bil_n ri_nj_diff_n" 
local nj_overall_c "ri_nj_eng_c ri_nj_bil_c ri_nj_diff_c" 

local va_ds_p16_n "ri_va_p16_eng_n ri_va_p16_bil_n ri_va_p16_diff_n"
local va_ds_p16_c "ri_va_p16_eng_c ri_va_p16_bil_c ri_va_p16_diff_c"

local va_ds_g16_n "ri_va_g16_eng_n ri_va_g16_bil_n ri_va_g16_diff_n"
local va_ds_g16_c "ri_va_g16_eng_c ri_va_g16_bil_c ri_va_g16_diff_c"

local nj_ds_p16_n "ri_nj_p16_eng_n ri_nj_p16_bil_n ri_nj_p16_diff_n"
local nj_ds_p16_c "ri_nj_p16_eng_c ri_nj_p16_bil_c ri_nj_p16_diff_c"

local nj_ds_g16_n "ri_nj_g16_eng_n ri_nj_g16_bil_n ri_nj_g16_diff_n"
local nj_ds_g16_c "ri_nj_g16_eng_c ri_nj_g16_bil_c ri_nj_g16_diff_c"
	
local va_group_n "ri_va_other_eng_n ri_va_other_bil_n ri_va_other_diff_n ri_va_mex_eng_n ri_va_mex_bil_n ri_va_mex_diff_n ri_va_mex_eng_h_n ri_va_mex_bil_h_n ri_va_mexother_diff_diff_n"
local va_group_c "ri_va_other_eng_c ri_va_other_bil_c ri_va_other_diff_c ri_va_mex_eng_c ri_va_mex_bil_c ri_va_mex_diff_c ri_va_mex_eng_h_c ri_va_mex_bil_h_c ri_va_mexother_diff_diff_c"

local nj_group_n "ri_nj_other_eng_n ri_nj_other_bil_n ri_nj_other_diff_n ri_nj_mex_eng_n ri_nj_mex_bil_n ri_nj_mex_diff_n ri_nj_pr_eng_n ri_nj_pr_bil_n ri_nj_pr_diff_n ri_nj_mex_eng_h_n ri_nj_mex_bil_h_n ri_nj_pr_eng_h_n ri_nj_pr_bil_h_n ri_nj_mexother_diff_diff_n ri_nj_prother_diff_diff_n ri_nj_mexpr_diff_diff_n"
local nj_group_c "ri_nj_other_eng_c ri_nj_other_bil_c ri_nj_other_diff_c ri_nj_mex_eng_c ri_nj_mex_bil_c ri_nj_mex_diff_c ri_nj_pr_eng_c ri_nj_pr_bil_c ri_nj_pr_diff_c ri_nj_mex_eng_h_c ri_nj_mex_bil_h_c ri_nj_pr_eng_h_c ri_nj_pr_bil_h_c ri_nj_mexother_diff_diff_c ri_nj_prother_diff_diff_c ri_nj_mexpr_diff_diff_c"

local nj_multi_n "ri_nj_shh_eng_n ri_nj_shh_bil_n ri_nj_shh_diff_n ri_nj_multi_eng_n ri_nj_multi_bil_n ri_nj_multi_diff_n ri_nj_multi_eng_h_n ri_nj_multi_bil_h_n ri_nj_multi_diff_diff_n" 
local nj_multi_c "ri_nj_shh_eng_c ri_nj_shh_bil_c ri_nj_shh_diff_c ri_nj_multi_eng_c ri_nj_multi_bil_c ri_nj_multi_diff_c ri_nj_multi_eng_h_c ri_nj_multi_bil_h_c ri_nj_multi_diff_diff_c" 

local nj_spill_n "ri_nj_direct_eng_n ri_nj_direct_bil_n ri_nj_direct_diff_n ri_nj_spill_eng_n ri_nj_spill_bil_n ri_nj_spill_diff_n" 
local nj_spill_c "ri_nj_direct_eng_c ri_nj_direct_bil_c ri_nj_direct_diff_c ri_nj_spill_eng_c ri_nj_spill_bil_c ri_nj_spill_diff_c"

tempname ri_estimates 
postfile `ri_estimates' ri_iter `va_overall_n' `va_overall_c' `nj_overall_n' `nj_overall_c' `va_ds_p16_n' `va_ds_p16_c' `va_ds_g16_n' `va_ds_g16_c' `nj_ds_p16_n' `nj_ds_p16_c' `nj_ds_g16_n' `nj_ds_g16_c' `va_group_n' `va_group_c' `nj_group_n' `nj_group_c' `nj_multi_n' `nj_multi_c' `nj_spill_n' `nj_spill_c' using ri_estimates.dta, replace		

* Start timer
timer on 1

*Set iterations
global ri_iterations = 1000

forval ri = 1/$ri_iterations {

*========================
di "RI ITERATION: `ri'"
*========================


*----------------------------------
* Randomization
*----------------------------------
use "Replication Data - NJ & VA 2015.dta", clear

gen assign_`ri' = .

	
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* VA Latino
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~

* Set the sort seed with random number generator.
set sortseed 339194

* Change seed by RI iteration 
local seed = 432909 + `ri'

* Block by household size & traditional rolloff vs consistent recent voters
local blockvars "g14 multi_hh"

* Balance by age, gender, vote history
local balancevars "electiondayage female g13 g12 g10 g08"


* Randomize ado - Wilks Lambda version (as of 10-09-2013)
	* Max tries = 10
	* acceptance threshold = p>0.8
randomize if hh_tag==1 & strata==1, groups(3) balance(`balancevars') block(`blockvars') seed(`seed') minruns(1) maxruns(10) jointp(0.8) replace

* Convert assignment variable in assigned groups 
bys mail_hhid: egen assign_tag_`ri'=max(_assignment)
replace assign_`ri' = 0 if assign_tag_`ri'==1 & assign_tag_`ri'!=.
replace assign_`ri' = 1 if assign_tag_`ri'==2 & assign_tag_`ri'!=.
replace assign_`ri' = 2 if assign_tag_`ri'==3 & assign_tag_`ri'!=.

tab assign_`ri' strata,  missing

drop _assignment assign_tag_`ri'


*~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* NJ Other Latino
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~

* Set the sort seed with random number generator.
set sortseed 307029

* Change seed by RI iteration 
local seed = 432909 + `ri'

* Block by household size & traditional rolloff vs consistent recent voters
local blockvars "g14"

* Balance by age, gender, vote history
local balancevars "electiondayage female g13 g12 g10 g08"

* Randomize ado - Wilks Lambda version (as of 10-09-2013)
	* Max tries = 10
	* acceptance threshold = p>0.8
randomize if hh_tag==1 & strata==3, groups(7) balance(`balancevars') block(`blockvars') seed(`seed') minruns(1) maxruns(10) jointp(0.8) replace

* Convert assignment variable in assigned groups 
replace assign_`ri' = 1 if _assignment==1 & strata==3
replace assign_`ri' = 2 if _assignment==2 & strata==3
replace assign_`ri' = 0 if _assignment>=3 & strata==3

tab assign_`ri' strata,  missing

drop _assignment



*~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* NJ Mexican
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~

* Set the sort seed with random number generator.
set sortseed 137667

* Change seed by RI iteration 
local seed = 432909 + `ri'

* Block by household size & traditional rolloff vs consistent recent voters
local blockvars "g14"

* Balance by age, gender, vote history
local balancevars "electiondayage female g13 g12 g10 g08"

* Randomize ado - Wilks Lambda version (as of 10-09-2013)
	* Max tries = 10
	* acceptance threshold = p>0.8
randomize if hh_tag==1 & strata==4, groups(6) balance(`balancevars') block(`blockvars') seed(`seed') minruns(1) maxruns(10) jointp(0.8) replace

* Convert assignment variable in assigned groups 
replace assign_`ri' = 1 if _assignment==1 & strata==4
replace assign_`ri' = 2 if _assignment==2 & strata==4
replace assign_`ri' = 0 if _assignment>=3 & strata==4

tab assign_`ri' strata,  missing

drop _assignment


*~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* NJ Puerto Rican
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~

* Set the sort seed with random number generator.
set sortseed 229125

* Change seed by RI iteration 
local seed = 432909 + `ri'

* Block by household size & traditional rolloff vs consistent recent voters
local blockvars "g14"

* Balance by age, gender, vote history
local balancevars "electiondayage female g13 g12 g10 g08"

* Randomize ado - Wilks Lambda version (as of 10-09-2013)
	* Max tries = 10
	* acceptance threshold = p>0.8
randomize if hh_tag==1 & strata==5, groups(7) balance(`balancevars') block(`blockvars') seed(`seed') minruns(1) maxruns(10) jointp(0.8) replace

* Convert assignment variable in assigned groups 
replace assign_`ri' = 1 if _assignment==1 & strata==5
replace assign_`ri' = 2 if _assignment==2 & strata==5
replace assign_`ri' = 0 if _assignment>=3 & strata==5

tab assign_`ri' strata,  missing

drop _assignment


*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Treatment Indicators for Latinos
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

* Create treatment dummies
gen eng_tx = assign_`ri' == 1 
gen bil_tx = assign_`ri' == 2 
	label var eng_tx "English Treatment"
	label var bil_tx "Bilingual Treatment"

* Interactions for Nation of Origin Strata
gen mex_eng = eng_tx * mexican
gen pr_eng = eng_tx * p_rican
gen mex_bil = bil_tx * mexican
gen pr_bil = bil_tx * p_rican

	label var mexican "Catalist: Mexican American"
	label var p_rican "Catalist: Puerto Rican"
	label var mex_eng "English * Mexican American"
	label var mex_bil "Bilingual * Mexican American"
	label var pr_eng "English * Puerto Rican"
	label var pr_bil "Bilingual * Puerto Rican"

* Interactions with multi target households
gen multi_eng = eng_tx * multi_hh
gen multi_bil = bil_tx * multi_hh
	label var multi_eng "eng_tx * Multi-target HH"
	label var multi_bil "bil_tx * Multi-target HH"

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Weights by Probability of Assignment
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Necessary because blocks do not have equal probability of assignment
* Weights are inverse probability of assignment

gen nj_wgt = .
replace nj_wgt = 1/(5/7) if inlist(strata,3,5) & assign_`ri'==0	
replace nj_wgt = 1/(1/7) if inlist(strata,3,5) & assign_`ri'==1	
replace nj_wgt = 1/(1/7) if inlist(strata,3,5) & assign_`ri'==2	

replace nj_wgt = 1/(4/6) if inlist(strata,4) & assign_`ri'==0	
replace nj_wgt = 1/(1/6) if inlist(strata,4) & assign_`ri'==1	
replace nj_wgt = 1/(1/6) if inlist(strata,4) & assign_`ri'==2	
	

*========================================
* ANALYSIS - SELECTED HH MEMBERS
*========================================
* Preserve randomization iteration for household spillover analysis
preserve

* Select relevant data
keep if hh_tag==1 | strata==1 /* All targets in HH mailed for VA Latino */

tab assign_`ri' strata, mi
 
*-------------------
* Virginia
*-------------------
* Set covariates & block variables for analysis
local balancevars "electiondayage female g14 g13 g12 g10 g08 multi_hh"

* eng_tx vs. bil_tx
reg g15 eng_tx bil_tx if strata==1 , robust cluster(mail_hhid)
	local ri_va_eng_n = _b[eng_tx]
	local ri_va_bil_n = _b[bil_tx]
	local ri_va_diff_n = _b[eng_tx]-_b[bil_tx]
reg g15 eng_tx bil_tx `balancevars' if strata==1 , robust cluster(mail_hhid)
	local ri_va_eng_c = _b[eng_tx]
	local ri_va_bil_c = _b[bil_tx]
	local ri_va_diff_c = _b[eng_tx]-_b[bil_tx]

/* Downstream in 2016 */
reg p16 eng_tx bil_tx if strata==1 , robust cluster(mail_hhid)
	local ri_va_p16_eng_n = _b[eng_tx]
	local ri_va_p16_bil_n = _b[bil_tx]
	local ri_va_p16_diff_n = _b[eng_tx]-_b[bil_tx]
reg p16 eng_tx bil_tx `balancevars' if strata==1 , robust cluster(mail_hhid)
	local ri_va_p16_eng_c = _b[eng_tx]
	local ri_va_p16_bil_c = _b[bil_tx]
	local ri_va_p16_diff_c = _b[eng_tx]-_b[bil_tx]

reg g16 eng_tx bil_tx if strata==1 , robust cluster(mail_hhid)
	local ri_va_g16_eng_n = _b[eng_tx]
	local ri_va_g16_bil_n = _b[bil_tx]
	local ri_va_g16_diff_n = _b[eng_tx]-_b[bil_tx]
reg g16 eng_tx bil_tx `balancevars' if strata==1 , robust cluster(mail_hhid)
	local ri_va_g16_eng_c = _b[eng_tx]
	local ri_va_g16_bil_c = _b[bil_tx]
	local ri_va_g16_diff_c = _b[eng_tx]-_b[bil_tx]


* Heterogeneity by nation of origin strata (without and with balance variables b/c blocked randomization)
reg g15 eng_tx mex_eng bil_tx mex_bil mexican if strata==1 , robust cluster(mail_hhid)
	local ri_va_other_eng_n = _b[eng_tx]
	local ri_va_other_bil_n = _b[bil_tx]
	local ri_va_other_diff_n = _b[eng_tx]-_b[bil_tx]
	local ri_va_mex_eng_n = _b[eng_tx] + _b[mex_eng]
	local ri_va_mex_bil_n = _b[bil_tx] + _b[mex_bil]
	local ri_va_mex_diff_n = (_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil])
	local ri_va_mex_eng_h_n = _b[mex_eng]
	local ri_va_mex_bil_h_n = _b[mex_bil]
	local ri_va_mexother_diff_diff_n = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))

reg g15 eng_tx mex_eng bil_tx mex_bil mexican `balancevars' if strata==1 , robust cluster(mail_hhid)
	local ri_va_other_eng_c = _b[eng_tx]
	local ri_va_other_bil_c = _b[bil_tx]
	local ri_va_other_diff_c = _b[eng_tx]-_b[bil_tx]
	local ri_va_mex_eng_c = _b[eng_tx] + _b[mex_eng]
	local ri_va_mex_bil_c = _b[bil_tx] + _b[mex_bil]
	local ri_va_mex_diff_c = (_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil])
	local ri_va_mex_eng_h_c = _b[mex_eng]
	local ri_va_mex_bil_h_c = _b[mex_bil]
	local ri_va_mexother_diff_diff_c = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))
		
*-------------------
* New Jersey
*-------------------
* Different probability of assignment across Other Latino, Mexican, Puerto Rican

* Set covariates & block variables for analysis (coding error in randomization omitted multi_hh)
local balancevars "electiondayage female g14 g13 g12 g10 g08"

* eng_tx vs. bil_tx
reg g15 eng_tx bil_tx mexican p_rican if strata>=2 [pweight=nj_wgt]
	local ri_nj_eng_n = _b[eng_tx]
	local ri_nj_bil_n = _b[bil_tx]
	local ri_nj_diff_n = _b[eng_tx]-_b[bil_tx]
reg g15 eng_tx bil_tx mexican p_rican `balancevars' if strata>=2 [pweight=nj_wgt]
	local ri_nj_eng_c = _b[eng_tx]
	local ri_nj_bil_c = _b[bil_tx]
	local ri_nj_diff_c = _b[eng_tx]-_b[bil_tx]

/* Downstream in 2016 */
reg p16 eng_tx bil_tx mexican p_rican if strata>=2 [pweight=nj_wgt]
	local ri_nj_p16_eng_n = _b[eng_tx]
	local ri_nj_p16_bil_n = _b[bil_tx]
	local ri_nj_p16_diff_n = _b[eng_tx]-_b[bil_tx]
reg p16 eng_tx bil_tx mexican p_rican `balancevars' if strata>=2 [pweight=nj_wgt]
	local ri_nj_p16_eng_c = _b[eng_tx]
	local ri_nj_p16_bil_c = _b[bil_tx]
	local ri_nj_p16_diff_c = _b[eng_tx]-_b[bil_tx]

reg g16 eng_tx bil_tx mexican p_rican if strata>=2 [pweight=nj_wgt]
	local ri_nj_g16_eng_n = _b[eng_tx]
	local ri_nj_g16_bil_n = _b[bil_tx]
	local ri_nj_g16_diff_n = _b[eng_tx]-_b[bil_tx]
reg g16 eng_tx bil_tx mexican p_rican `balancevars' if strata>=2 [pweight=nj_wgt]
	local ri_nj_g16_eng_c = _b[eng_tx]
	local ri_nj_g16_bil_c = _b[bil_tx]
	local ri_nj_g16_diff_c = _b[eng_tx]-_b[bil_tx]

* Heterogeneity by nation of origin strata (without and with balance variables b/c blocked randomization)
reg g15 eng_tx mex_eng pr_eng bil_tx mex_bil pr_bil mexican p_rican if strata>=2 [pweight=nj_wgt]
	local ri_nj_other_eng_n = _b[eng_tx]
	local ri_nj_other_bil_n = _b[bil_tx]
	local ri_nj_other_diff_n = _b[eng_tx]-_b[bil_tx]
	local ri_nj_mex_eng_n = _b[eng_tx] + _b[mex_eng]
	local ri_nj_mex_bil_n = _b[bil_tx] + _b[mex_bil]
	local ri_nj_mex_diff_n = (_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil])
	local ri_nj_pr_eng_n = _b[eng_tx] + _b[pr_eng]
	local ri_nj_pr_bil_n = _b[bil_tx] + _b[pr_bil]
	local ri_nj_pr_diff_n = (_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil])
	local ri_nj_mex_eng_h_n = _b[mex_eng]
	local ri_nj_mex_bil_h_n = _b[mex_bil]
	local ri_nj_pr_eng_h_n = _b[pr_eng]
	local ri_nj_pr_bil_h_n = _b[pr_bil]
	local ri_nj_mexother_diff_diff_n = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))
	local ri_nj_prother_diff_diff_n = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil]))
	local ri_nj_mexpr_diff_diff_n = ((_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil]))-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))
reg g15 eng_tx mex_eng pr_eng bil_tx mex_bil pr_bil mexican p_rican `balancevars' if strata>=2 [pweight=nj_wgt]
	local ri_nj_other_eng_c = _b[eng_tx]
	local ri_nj_other_bil_c = _b[bil_tx]
	local ri_nj_other_diff_c = _b[eng_tx]-_b[bil_tx]
	local ri_nj_mex_eng_c = _b[eng_tx] + _b[mex_eng]
	local ri_nj_mex_bil_c = _b[bil_tx] + _b[mex_bil]
	local ri_nj_mex_diff_c = (_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil])
	local ri_nj_pr_eng_c = _b[eng_tx] + _b[pr_eng]
	local ri_nj_pr_bil_c = _b[bil_tx] + _b[pr_bil]
	local ri_nj_pr_diff_c = (_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil])
	local ri_nj_mex_eng_h_c = _b[mex_eng]
	local ri_nj_mex_bil_h_c = _b[mex_bil]
	local ri_nj_pr_eng_h_c = _b[pr_eng]
	local ri_nj_pr_bil_h_c = _b[pr_bil]	
	local ri_nj_mexother_diff_diff_c = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))
	local ri_nj_prother_diff_diff_c = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil]))
	local ri_nj_mexpr_diff_diff_c = ((_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil]))-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))
	
* Heterogeneity multi-target households
reg g15 eng_tx multi_eng bil_tx multi_bil multi_hh mexican p_rican hhsize_3 hhsize_4 hhsize_5 if strata>=2 [pweight=nj_wgt]
	local ri_nj_shh_eng_n = _b[eng_tx]
	local ri_nj_shh_bil_n = _b[bil_tx]
	local ri_nj_shh_diff_n = _b[eng_tx]-_b[bil_tx]
	local ri_nj_multi_eng_n = _b[eng_tx] + _b[multi_eng]
	local ri_nj_multi_bil_n = _b[bil_tx] + _b[multi_bil]
	local ri_nj_multi_diff_n = (_b[eng_tx] + _b[multi_eng]) - (_b[bil_tx] + _b[multi_bil])
	local ri_nj_multi_eng_h_n = _b[multi_eng]
	local ri_nj_multi_bil_h_n = _b[multi_bil]
	local ri_nj_multi_diff_diff_n = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[multi_eng]) - (_b[bil_tx] + _b[multi_bil]))
reg g15 eng_tx multi_eng bil_tx multi_bil multi_hh mexican p_rican `balancevars' hhsize_3 hhsize_4 hhsize_5 if strata>=2 [pweight=nj_wgt]
	local ri_nj_shh_eng_c = _b[eng_tx]
	local ri_nj_shh_bil_c = _b[bil_tx]
	local ri_nj_shh_diff_c = _b[eng_tx]-_b[bil_tx]
	local ri_nj_multi_eng_c = _b[eng_tx] + _b[multi_eng]
	local ri_nj_multi_bil_c = _b[bil_tx] + _b[multi_bil]
	local ri_nj_multi_diff_c = (_b[eng_tx] + _b[multi_eng]) - (_b[bil_tx] + _b[multi_bil])
	local ri_nj_multi_eng_h_c = _b[multi_eng]
	local ri_nj_multi_bil_h_c = _b[multi_bil]
	local ri_nj_multi_diff_diff_c = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[multi_eng]) - (_b[bil_tx] + _b[multi_bil]))
	
	
*====================================================
* ANALYSIS - Direct & indirect ATE in multiperson HH
*====================================================

* Restore randomization iteration (and clear memory)
restore

* Select relevant data
keep if hh_count>=2 
drop if strata==1 /* All targets in HH mailed for VA Latino */


*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Treatment Indicators for Latinos
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

* Spillover treatments
bys mail_hhid: egen assign_hh_`ri' = max(assign_`ri')
tab assign_hh_`ri' strata
gen eng_tx_hh = assign_hh_`ri' == 1 
gen bil_tx_hh = assign_hh_`ri' == 2 



*-------------------
* New Jersey
*-------------------
* Different probability of assignment across Other Latino, Mexican, Puerto Rican
	local balancevars "electiondayage female g14 g13 g12 g10 g08"

* Direct Treatment - eng_tx vs. bil_tx
reg g15 eng_tx bil_tx mexican p_rican hhsize_3 hhsize_4 hhsize_5 if hh_tag==1 [pweight=nj_wgt]
	local ri_nj_direct_eng_n = _b[eng_tx]
	local ri_nj_direct_bil_n = _b[bil_tx]
	local ri_nj_direct_diff_n = _b[eng_tx]-_b[bil_tx]
reg g15 eng_tx bil_tx mexican p_rican hhsize_3 hhsize_4 hhsize_5 `balancevars' if hh_tag==1 [pweight=nj_wgt]
	local ri_nj_direct_eng_c = _b[eng_tx]
	local ri_nj_direct_bil_c = _b[bil_tx]
	local ri_nj_direct_diff_c = _b[eng_tx]-_b[bil_tx]

* Spillover - eng_tx vs. bil_tx
reg g15 eng_tx_hh bil_tx_hh mexican p_rican hhsize_3 hhsize_4 hhsize_5 if hh_tag==0 [pweight=nj_wgt], robust cluster(mail_hhid) 
	local ri_nj_spill_eng_n = _b[eng_tx_hh]
	local ri_nj_spill_bil_n = _b[bil_tx_hh]
	local ri_nj_spill_diff_n = _b[eng_tx_hh]-_b[bil_tx_hh]
reg g15 eng_tx_hh bil_tx_hh mexican p_rican hhsize_3 hhsize_4 hhsize_5 `balancevars' if hh_tag==0 [pweight=nj_wgt], robust cluster(mail_hhid) 
	local ri_nj_spill_eng_c = _b[eng_tx_hh]
	local ri_nj_spill_bil_c = _b[bil_tx_hh]
	local ri_nj_spill_diff_c = _b[eng_tx_hh]-_b[bil_tx_hh]

post `ri_estimates' (`ri') ///
					(`ri_va_eng_n') (`ri_va_bil_n') (`ri_va_diff_n') ///
					(`ri_va_eng_c') (`ri_va_bil_c') (`ri_va_diff_c') ///
					(`ri_nj_eng_n') (`ri_nj_bil_n') (`ri_nj_diff_n') ///
					(`ri_nj_eng_c') (`ri_nj_bil_c') (`ri_nj_diff_c') ///
					(`ri_va_p16_eng_n') (`ri_va_p16_bil_n') (`ri_va_p16_diff_n') ///
					(`ri_va_p16_eng_c') (`ri_va_p16_bil_c') (`ri_va_p16_diff_c') ///
					(`ri_va_g16_eng_n') (`ri_va_g16_bil_n') (`ri_va_g16_diff_n') ///
					(`ri_va_g16_eng_c') (`ri_va_g16_bil_c') (`ri_va_g16_diff_c') ///
					(`ri_nj_p16_eng_n') (`ri_nj_p16_bil_n') (`ri_nj_p16_diff_n') ///
					(`ri_nj_p16_eng_c') (`ri_nj_p16_bil_c') (`ri_nj_p16_diff_c') ///
					(`ri_nj_g16_eng_n') (`ri_nj_g16_bil_n') (`ri_nj_g16_diff_n') ///
					(`ri_nj_g16_eng_c') (`ri_nj_g16_bil_c') (`ri_nj_g16_diff_c') ///
					(`ri_va_other_eng_n') (`ri_va_other_bil_n') (`ri_va_other_diff_n') ///
					(`ri_va_mex_eng_n') (`ri_va_mex_bil_n') (`ri_va_mex_diff_n') ///
					(`ri_va_mex_eng_h_n') (`ri_va_mex_bil_h_n') (`ri_va_mexother_diff_diff_n') ///
					(`ri_va_other_eng_c') (`ri_va_other_bil_c') (`ri_va_other_diff_c') ///
					(`ri_va_mex_eng_c') (`ri_va_mex_bil_c') (`ri_va_mex_diff_c') ///
					(`ri_va_mex_eng_h_c') (`ri_va_mex_bil_h_c') (`ri_va_mexother_diff_diff_c') ///
					(`ri_nj_other_eng_n') (`ri_nj_other_bil_n') (`ri_nj_other_diff_n') ///
					(`ri_nj_mex_eng_n') (`ri_nj_mex_bil_n') (`ri_nj_mex_diff_n') ///
					(`ri_nj_pr_eng_n') (`ri_nj_pr_bil_n') (`ri_nj_pr_diff_n') ///
					(`ri_nj_mex_eng_h_n') (`ri_nj_mex_bil_h_n') ///
					(`ri_nj_pr_eng_h_n') (`ri_nj_pr_bil_h_n') ///
					(`ri_nj_mexother_diff_diff_n') (`ri_nj_prother_diff_diff_n') (`ri_nj_mexpr_diff_diff_n') ///
					(`ri_nj_other_eng_c') (`ri_nj_other_bil_c') (`ri_nj_other_diff_c') ///
					(`ri_nj_mex_eng_c') (`ri_nj_mex_bil_c') (`ri_nj_mex_diff_c') ///
					(`ri_nj_pr_eng_c') (`ri_nj_pr_bil_c') (`ri_nj_pr_diff_c') ///
					(`ri_nj_mex_eng_h_c') (`ri_nj_mex_bil_h_c') ///
					(`ri_nj_pr_eng_h_c') (`ri_nj_pr_bil_h_c') ///
					(`ri_nj_mexother_diff_diff_c') (`ri_nj_prother_diff_diff_c') (`ri_nj_mexpr_diff_diff_c') ///
					(`ri_nj_shh_eng_n') (`ri_nj_shh_bil_n') (`ri_nj_shh_diff_n') ///
					(`ri_nj_multi_eng_n') (`ri_nj_multi_bil_n') (`ri_nj_multi_diff_n') ///
					(`ri_nj_multi_eng_h_n') (`ri_nj_multi_bil_h_n') (`ri_nj_multi_diff_diff_n') ///
					(`ri_nj_shh_eng_c') (`ri_nj_shh_bil_c') (`ri_nj_shh_diff_c') ///
					(`ri_nj_multi_eng_c') (`ri_nj_multi_bil_c') (`ri_nj_multi_diff_c') ///
					(`ri_nj_multi_eng_h_c') (`ri_nj_multi_bil_h_c') (`ri_nj_multi_diff_diff_c') ///
					(`ri_nj_direct_eng_n') (`ri_nj_direct_bil_n') (`ri_nj_direct_diff_n') ///
					(`ri_nj_spill_eng_n') (`ri_nj_spill_bil_n') (`ri_nj_spill_diff_n') ///
					(`ri_nj_direct_eng_c') (`ri_nj_direct_bil_c') (`ri_nj_direct_diff_c') ///
					(`ri_nj_spill_eng_c') (`ri_nj_spill_bil_c') (`ri_nj_spill_diff_c')

	
}
postclose `ri_estimates'
	
* Check time elapsed
timer off 1
timer list 1
disp "Elapsed time (in hours) = " r(t1)/(60*60) 	
	
cd "<filepath>"
	
*----------------------------------------
* Capture Observed ATEs
*----------------------------------------
use "Replication Data - NJ & VA 2015.dta", clear

drop if hispanic!=1

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Treatment Indicators for Latinos
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

gen eng_tx = assign_tx == 1 
gen bil_tx = assign_tx == 2 
	label var eng_tx "English Treatment"
	label var bil_tx "Bilingual Treatment"


* Interactions for Nation of Origin Strata
replace mexican = ethnicity == "Mexican" if state=="VA"
gen mex_eng = eng_tx * mexican
gen pr_eng = eng_tx * p_rican
gen mex_bil = bil_tx * mexican
gen pr_bil = bil_tx * p_rican

	label var mexican "Catalist: Mexican American"
	label var p_rican "Catalist: Puerto Rican"
	label var mex_eng "English * Mexican American"
	label var mex_bil "Bilingual * Mexican American"
	label var pr_eng "English * Puerto Rican"
	label var pr_bil "Bilingual * Puerto Rican"

	
* Interactions with multi target households
gen multi_eng = eng_tx * multi_hh
gen multi_bil = bil_tx * multi_hh
	label var multi_eng "eng_tx * Multi-target HH"
	label var multi_bil "bil_tx * Multi-target HH"

* Covariates for Household size
/* omit HH=2 as base -- gen hhsize_2 = hh_count== 2 */
gen hhsize_3 = hh_count== 3
gen hhsize_4 = hh_count== 4
gen hhsize_5 = hh_count>= 5
	label var hhsize_3 "HH of 3 people"
	label var hhsize_4 "HH of 4 people"
	label var hhsize_5 "HH 5 or more people"


*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Weights by Probability of Assignment
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Necessary because blocks do not have equal probability of assignment
* Weights are inverse probability of assignment

gen nj_wgt = .
replace nj_wgt = 1/(5/7) if inlist(strata,3,5) & assign_`ri'==0	
replace nj_wgt = 1/(1/7) if inlist(strata,3,5) & assign_`ri'==1	
replace nj_wgt = 1/(1/7) if inlist(strata,3,5) & assign_`ri'==2	

replace nj_wgt = 1/(4/6) if inlist(strata,4) & assign_`ri'==0	
replace nj_wgt = 1/(1/6) if inlist(strata,4) & assign_`ri'==1	
replace nj_wgt = 1/(1/6) if inlist(strata,4) & assign_`ri'==2	
	

*========================================
* ANALYSIS - SELECTED HH MEMBERS
*========================================
* Preserve randomization iteration for household spillover analysis
preserve

* Select relevant data
keep if hh_tag==1 | strata==1 /* All targets in HH mailed for VA Latino */

 
*-------------------
* Virginia
*-------------------
* Set covariates & block variables for analysis
local balancevars "electiondayage female g14 g13 g12 g10 g08 multi_hh"

* eng_tx vs. bil_tx
reg g15 eng_tx bil_tx if strata==1 , robust cluster(mail_hhid)
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_va_eng_n = _b[eng_tx]
	local obs_va_bil_n = _b[bil_tx]
	local obs_va_diff_n = _b[eng_tx]-_b[bil_tx]
reg g15 eng_tx bil_tx `balancevars' if strata==1 , robust cluster(mail_hhid)
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_va_eng_c = _b[eng_tx]
	local obs_va_bil_c = _b[bil_tx]
	local obs_va_diff_c = _b[eng_tx]-_b[bil_tx]

/* Downstream in 2016 */
reg p16 eng_tx bil_tx if strata==1 , robust cluster(mail_hhid)
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_va_p16_eng_n = _b[eng_tx]
	local obs_va_p16_bil_n = _b[bil_tx]
	local obs_va_p16_diff_n = _b[eng_tx]-_b[bil_tx]

reg p16 eng_tx bil_tx `balancevars' if strata==1 , robust cluster(mail_hhid)
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_va_p16_eng_c = _b[eng_tx]
	local obs_va_p16_bil_c = _b[bil_tx]
	local obs_va_p16_diff_c = _b[eng_tx]-_b[bil_tx]

reg g16 eng_tx bil_tx if strata==1 , robust cluster(mail_hhid)
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_va_g16_eng_n = _b[eng_tx]
	local obs_va_g16_bil_n = _b[bil_tx]
	local obs_va_g16_diff_n = _b[eng_tx]-_b[bil_tx]

reg g16 eng_tx bil_tx `balancevars' if strata==1 , robust cluster(mail_hhid)
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_va_g16_eng_c = _b[eng_tx]
	local obs_va_g16_bil_c = _b[bil_tx]
	local obs_va_g16_diff_c = _b[eng_tx]-_b[bil_tx]


* Heterogeneity by nation of origin strata (without and with balance variables b/c blocked randomization)
reg g15 eng_tx mex_eng bil_tx mex_bil mexican if strata==1 , robust cluster(mail_hhid)
	* CATE
	lincom eng_tx + mex_eng
	lincom bil_tx + mex_bil
	* Diff in ATE
	lincom eng_tx - bil_tx
	lincom (eng_tx + mex_eng) - (bil_tx + mex_bil)
	lincom (eng_tx - bil_tx) - ((eng_tx + mex_eng) - (bil_tx + mex_bil))

	local obs_va_other_eng_n = _b[eng_tx]
	local obs_va_other_bil_n = _b[bil_tx]
	local obs_va_other_diff_n = _b[eng_tx]-_b[bil_tx]
	local obs_va_mex_eng_n = _b[eng_tx] + _b[mex_eng]
	local obs_va_mex_bil_n = _b[bil_tx] + _b[mex_bil]
	local obs_va_mex_diff_n = (_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil])
	local obs_va_mex_eng_h_n = _b[mex_eng]
	local obs_va_mex_bil_h_n = _b[mex_bil]
	local obs_va_mexother_diff_diff_n = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))

reg g15 eng_tx mex_eng bil_tx mex_bil mexican `balancevars' if strata==1 , robust cluster(mail_hhid)
	* CATE
	lincom eng_tx + mex_eng
	lincom bil_tx + mex_bil
	* Diff in ATE
	lincom eng_tx - bil_tx
	lincom (eng_tx + mex_eng) - (bil_tx + mex_bil)
	lincom (eng_tx - bil_tx) - ((eng_tx + mex_eng) - (bil_tx + mex_bil))

	local obs_va_other_eng_c = _b[eng_tx]
	local obs_va_other_bil_c = _b[bil_tx]
	local obs_va_other_diff_c = _b[eng_tx]-_b[bil_tx]
	local obs_va_mex_eng_c = _b[eng_tx] + _b[mex_eng]
	local obs_va_mex_bil_c = _b[bil_tx] + _b[mex_bil]
	local obs_va_mex_diff_c = (_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil])
	local obs_va_mex_eng_h_c = _b[mex_eng]
	local obs_va_mex_bil_h_c = _b[mex_bil]
	local obs_va_mexother_diff_diff_c = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))
		
*-------------------
* New Jersey
*-------------------
* Different probability of assignment across Other Latino, Mexican, Puerto Rican

* Set covariates & block variables for analysis (coding error in randomization omitted multi_hh)
local balancevars "electiondayage female g14 g13 g12 g10 g08"

* eng_tx vs. bil_tx
reg g15 eng_tx bil_tx mexican p_rican if strata>=2 [pweight=nj_wgt]
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_eng_n = _b[eng_tx]
	local obs_nj_bil_n = _b[bil_tx]
	local obs_nj_diff_n = _b[eng_tx]-_b[bil_tx]
reg g15 eng_tx bil_tx mexican p_rican `balancevars' if strata>=2 [pweight=nj_wgt]
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_eng_c = _b[eng_tx]
	local obs_nj_bil_c = _b[bil_tx]
	local obs_nj_diff_c = _b[eng_tx]-_b[bil_tx]

/* Downstream in 2016 */
reg p16 eng_tx bil_tx mexican p_rican if strata>=2 [pweight=nj_wgt]
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_p16_eng_n = _b[eng_tx]
	local obs_nj_p16_bil_n = _b[bil_tx]
	local obs_nj_p16_diff_n = _b[eng_tx]-_b[bil_tx]
reg p16 eng_tx bil_tx mexican p_rican `balancevars' if strata>=2 [pweight=nj_wgt]
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_p16_eng_c = _b[eng_tx]
	local obs_nj_p16_bil_c = _b[bil_tx]
	local obs_nj_p16_diff_c = _b[eng_tx]-_b[bil_tx]

reg g16 eng_tx bil_tx mexican p_rican if strata>=2 [pweight=nj_wgt]
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_g16_eng_n = _b[eng_tx]
	local obs_nj_g16_bil_n = _b[bil_tx]
	local obs_nj_g16_diff_n = _b[eng_tx]-_b[bil_tx]
reg g16 eng_tx bil_tx mexican p_rican `balancevars' if strata>=2 [pweight=nj_wgt]
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_g16_eng_c = _b[eng_tx]
	local obs_nj_g16_bil_c = _b[bil_tx]
	local obs_nj_g16_diff_c = _b[eng_tx]-_b[bil_tx]


* Heterogeneity by nation of origin strata (without and with balance variables b/c blocked randomization)
reg g15 eng_tx mex_eng pr_eng bil_tx mex_bil pr_bil mexican p_rican if strata>=2 [pweight=nj_wgt]
	* CATE
	lincom eng_tx + mex_eng
	lincom eng_tx + pr_eng
	lincom bil_tx + mex_bil
	lincom bil_tx + pr_bil
	* Diff in ATE
	lincom eng_tx - bil_tx
	lincom (eng_tx + mex_eng) - (bil_tx + mex_bil)
	lincom (eng_tx + pr_eng) - (bil_tx + pr_bil)
	lincom (eng_tx - bil_tx) - ((eng_tx + mex_eng) - (bil_tx + mex_bil)) 
	lincom (eng_tx - bil_tx) - ((eng_tx + pr_eng) - (bil_tx + pr_bil))
	lincom ((eng_tx + pr_eng) - (bil_tx + pr_bil)) - ((eng_tx + mex_eng) - (bil_tx + mex_bil)) 

	local obs_nj_other_eng_n = _b[eng_tx]
	local obs_nj_other_bil_n = _b[bil_tx]
	local obs_nj_other_diff_n = _b[eng_tx]-_b[bil_tx]
	local obs_nj_mex_eng_n = _b[eng_tx] + _b[mex_eng]
	local obs_nj_mex_bil_n = _b[bil_tx] + _b[mex_bil]
	local obs_nj_mex_diff_n = (_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil])
	local obs_nj_pr_eng_n = _b[eng_tx] + _b[pr_eng]
	local obs_nj_pr_bil_n = _b[bil_tx] + _b[pr_bil]
	local obs_nj_pr_diff_n = (_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil])
	local obs_nj_mex_eng_h_n = _b[mex_eng]
	local obs_nj_mex_bil_h_n = _b[mex_bil]
	local obs_nj_pr_eng_h_n = _b[pr_eng]
	local obs_nj_pr_bil_h_n = _b[pr_bil]
	local obs_nj_mexother_diff_diff_n = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))
	local obs_nj_prother_diff_diff_n = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil]))
	local obs_nj_mexpr_diff_diff_n = ((_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil]))-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))

reg g15 eng_tx mex_eng pr_eng bil_tx mex_bil pr_bil mexican p_rican `balancevars' if strata>=2 [pweight=nj_wgt]
	* CATE
	lincom eng_tx + mex_eng
	lincom eng_tx + pr_eng
	lincom bil_tx + mex_bil
	lincom bil_tx + pr_bil
	* Diff in ATE
	lincom eng_tx - bil_tx
	lincom (eng_tx + mex_eng) - (bil_tx + mex_bil)
	lincom (eng_tx + pr_eng) - (bil_tx + pr_bil)
	lincom (eng_tx - bil_tx) - ((eng_tx + mex_eng) - (bil_tx + mex_bil)) 
	lincom (eng_tx - bil_tx) - ((eng_tx + pr_eng) - (bil_tx + pr_bil))
	lincom ((eng_tx + pr_eng) - (bil_tx + pr_bil)) - ((eng_tx + mex_eng) - (bil_tx + mex_bil)) 

	local obs_nj_other_eng_c = _b[eng_tx]
	local obs_nj_other_bil_c = _b[bil_tx]
	local obs_nj_other_diff_c = _b[eng_tx]-_b[bil_tx]
	local obs_nj_mex_eng_c = _b[eng_tx] + _b[mex_eng]
	local obs_nj_mex_bil_c = _b[bil_tx] + _b[mex_bil]
	local obs_nj_mex_diff_c = (_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil])
	local obs_nj_pr_eng_c = _b[eng_tx] + _b[pr_eng]
	local obs_nj_pr_bil_c = _b[bil_tx] + _b[pr_bil]
	local obs_nj_pr_diff_c = (_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil])
	local obs_nj_mex_eng_h_c = _b[mex_eng]
	local obs_nj_mex_bil_h_c = _b[mex_bil]
	local obs_nj_pr_eng_h_c = _b[pr_eng]
	local obs_nj_pr_bil_h_c = _b[pr_bil]	
	local obs_nj_mexother_diff_diff_c = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))
	local obs_nj_prother_diff_diff_c = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil]))
	local obs_nj_mexpr_diff_diff_c = ((_b[eng_tx] + _b[pr_eng]) - (_b[bil_tx] + _b[pr_bil]))-((_b[eng_tx] + _b[mex_eng]) - (_b[bil_tx] + _b[mex_bil]))
	
* Heterogeneity multi-target households
reg g15 eng_tx multi_eng bil_tx multi_bil multi_hh mexican p_rican hhsize_3 hhsize_4 hhsize_5 if strata>=2 [pweight=nj_wgt]
	* CATE
	lincom eng_tx + multi_eng
	lincom bil_tx + multi_bil
	* Diff in ATE
	lincom eng_tx - bil_tx
	lincom (eng_tx + multi_eng) - (bil_tx + multi_bil)
	lincom (eng_tx - bil_tx) - ((eng_tx + multi_eng) - (bil_tx + multi_bil))

	local obs_nj_shh_eng_n = _b[eng_tx]
	local obs_nj_shh_bil_n = _b[bil_tx]
	local obs_nj_shh_diff_n = _b[eng_tx]-_b[bil_tx]
	local obs_nj_multi_eng_n = _b[eng_tx] + _b[multi_eng]
	local obs_nj_multi_bil_n = _b[bil_tx] + _b[multi_bil]
	local obs_nj_multi_diff_n = (_b[eng_tx] + _b[multi_eng]) - (_b[bil_tx] + _b[multi_bil])
	local obs_nj_multi_eng_h_n = _b[multi_eng]
	local obs_nj_multi_bil_h_n = _b[multi_bil]
	local obs_nj_multi_diff_diff_n = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[multi_eng]) - (_b[bil_tx] + _b[multi_bil]))
reg g15 eng_tx multi_eng bil_tx multi_bil multi_hh mexican p_rican `balancevars' hhsize_3 hhsize_4 hhsize_5 if strata>=2 [pweight=nj_wgt]
	* CATE
	lincom eng_tx + multi_eng
	lincom bil_tx + multi_bil
	* Diff in ATE
	lincom eng_tx - bil_tx
	lincom (eng_tx + multi_eng) - (bil_tx + multi_bil)
	lincom (eng_tx - bil_tx) - ((eng_tx + multi_eng) - (bil_tx + multi_bil))

	
	local obs_nj_shh_eng_c = _b[eng_tx]
	local obs_nj_shh_bil_c = _b[bil_tx]
	local obs_nj_shh_diff_c = _b[eng_tx]-_b[bil_tx]
	local obs_nj_multi_eng_c = _b[eng_tx] + _b[multi_eng]
	local obs_nj_multi_bil_c = _b[bil_tx] + _b[multi_bil]
	local obs_nj_multi_diff_c = (_b[eng_tx] + _b[multi_eng]) - (_b[bil_tx] + _b[multi_bil])
	local obs_nj_multi_eng_h_c = _b[multi_eng]
	local obs_nj_multi_bil_h_c = _b[multi_bil]
	local obs_nj_multi_diff_diff_c = (_b[eng_tx]-_b[bil_tx])-((_b[eng_tx] + _b[multi_eng]) - (_b[bil_tx] + _b[multi_bil]))
	
	
*====================================================
* ANALYSIS - Direct & indirect ATE in multiperson HH
*====================================================

* Restore randomization iteration (and clear memory)
restore

* Select relevant data
keep if hh_count>=2 
drop if strata==1 /* All targets in HH mailed for VA Latino */

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Treatment Indicators for Latinos
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


* Direct Treatments
*gen eng_tx = assign_tx == 1 
*gen bil_tx = assign_tx == 2 


* Spillover treatments
bys mail_hhid: egen assign_tx_hh = max(assign_tx)
tab assign_tx_hh strata
gen eng_tx_hh = assign_tx_hh == 1 
gen bil_tx_hh = assign_tx_hh == 2 

* Covariates for Household size
/* omit HH=2 as base -- gen hhsize_2 = hh_count== 2 */
*gen hhsize_3 = hh_count== 3
*gen hhsize_4 = hh_count== 4
*gen hhsize_5 = hh_count>= 5



*-------------------
* New Jersey
*-------------------
* Different probability of assignment across Other Latino, Mexican, Puerto Rican
	local balancevars "electiondayage female g14 g13 g12 g10 g08"

* Direct Treatment - eng_tx vs. bil_tx
reg g15 eng_tx bil_tx mexican p_rican hhsize_3 hhsize_4 hhsize_5 if hh_tag==1 [pweight=nj_wgt]
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_direct_eng_n = _b[eng_tx]
	local obs_nj_direct_bil_n = _b[bil_tx]
	local obs_nj_direct_diff_n = _b[eng_tx]-_b[bil_tx]
reg g15 eng_tx bil_tx mexican p_rican hhsize_3 hhsize_4 hhsize_5 `balancevars' if hh_tag==1 [pweight=nj_wgt]
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_direct_eng_c = _b[eng_tx]
	local obs_nj_direct_bil_c = _b[bil_tx]
	local obs_nj_direct_diff_c = _b[eng_tx]-_b[bil_tx]

* Spillover - eng_tx vs. bil_tx
reg g15 eng_tx_hh bil_tx_hh mexican p_rican hhsize_3 hhsize_4 hhsize_5 if hh_tag==0 [pweight=nj_wgt], robust cluster(mail_hhid) 
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_spill_eng_n = _b[eng_tx_hh]
	local obs_nj_spill_bil_n = _b[bil_tx_hh]
	local obs_nj_spill_diff_n = _b[eng_tx_hh]-_b[bil_tx_hh]
reg g15 eng_tx_hh bil_tx_hh mexican p_rican hhsize_3 hhsize_4 hhsize_5 `balancevars' if hh_tag==0 [pweight=nj_wgt], robust cluster(mail_hhid) 
	* Diff in ATE
	lincom eng_tx - bil_tx

	local obs_nj_spill_eng_c = _b[eng_tx_hh]
	local obs_nj_spill_bil_c = _b[bil_tx_hh]
	local obs_nj_spill_diff_c = _b[eng_tx_hh]-_b[bil_tx_hh]


	
*-----------------------------------
* Calculate RI p-values 
*-----------------------------------
cd "<filepath>"

* Open file of RI estimates
use ri_estimates.dta, clear

* quantity is local that is list of all quantities of interest
local quantity "va_eng va_bil va_diff va_p16_eng va_p16_bil va_p16_diff va_g16_eng va_g16_bil va_g16_diff va_other_eng va_other_bil va_other_diff va_mex_eng va_mex_bil va_mex_diff va_mex_eng_h va_mex_bil_h va_mexother_diff_diff nj_eng nj_bil nj_diff nj_p16_eng nj_p16_bil nj_p16_diff nj_g16_eng nj_g16_bil nj_g16_diff nj_other_eng nj_other_bil nj_other_diff nj_mex_eng nj_mex_bil nj_mex_diff nj_pr_eng nj_pr_bil nj_pr_diff nj_mex_eng_h nj_mex_bil_h nj_pr_eng_h nj_pr_bil_h nj_mexother_diff_diff nj_prother_diff_diff nj_mexpr_diff_diff nj_shh_eng nj_shh_bil nj_shh_diff nj_multi_eng nj_multi_bil nj_multi_diff nj_multi_eng_h nj_multi_bil_h nj_multi_diff_diff nj_direct_eng nj_direct_bil nj_direct_diff nj_spill_eng nj_spill_bil nj_spill_diff"

* RI p-values from model without and with covariates (two-sided hypothesis test)
foreach model in "n" "c" {
  foreach outcome in `quantity' {
	gen c_`outcome'_`model' = abs(ri_`outcome'_`model')> abs(`obs_`outcome'_`model'') 
	label var c_`outcome'_`model' "Cases w/ `outcome'_`model' > observed value (two-sided)"
	qui sum c_`outcome'_`model'
		local p_`outcome'_`model' = r(mean)
	display "estimate for `outcome'_`model' = `obs_`outcome'_`model''"
	display "`outcome'_`model' p-value =  `p_`outcome'_`model''"
		
	di "-------------------------------------"
	
	
  }
}


log close
